Wave nucleation rate in excitable systems in the low noise limit 
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Motivated by recent experiments on intracellular calcium dynamics, we study the general issue of 
fluctuation-induced nucleation of waves in excitable media. We utilize a stochastic Fitzhugh-Nagumo 
model for this study, a spatially-extended non-potential pair of equations driven by thermal (i.e. 
white) noise. The nucleation rate is determined by finding the most probable escape path via 
minimization of an action related to the deviation of the fields from their deterministic trajectories. 
Our results pave the way both for studies of more realistic models of calcium dynamics as well as 
of nucleation phenomena in other non-equilibrium pattern-forming processes. 



One very important class of non-equilibrium spatially- 
extended systems is that of excitable media. In these, a 
quiescent state is linearly stable but nonlinear waves can 
nonetheless propagate without decaying. These waves 
can be generated by above-threshold local perturbations 
and they also can become self-sustaining in the form of 
rotating spirals0. Examples of excitable media include 
many biological systems such as the cAMP waves seen in 
Dictyostelium amoebae aggregation 0, electrical waves 
in cardiac and neural tissue[3j and, primary for our focus 
here, intracellular calcium waves 

Most excitable systems are sufficiently macroscopic as 
to render unimportant the role of thermodynamic fluc- 
tuations and allow for a description in terms of deter- 
ministic pde models. For these cases, noise effects can 
still be studied via the imposition of external variation 
in time and/or space (e.g. by varying illumination in a 
light-sensitive BZ reaction ja, 0|); however, there is no 
need to include noise in a description of the "natural" 
version of these systems. This is manifestly not the case 




FIG. 1: Left: Phase space diagram of the Fitzhugh-Nagumo 
system without noise and without the diffusive term. The 
stable equilibrium point is O. The f(u, v) — 3u — v? — v 
nullcline is the S shaped solid line with minimum at M and 
the g = u + (1 + 7) nullcline is the thin solid line. A small 
perturbation of O leads to a large excursion in phase space. 
The dashed line shows the phase-space excursion of a point 
during the propagation of a single wave. Parameter values are 
7 = .5 and e = 0.05. Right: Typical propagating wave in an 
excitable medium. Solid line: u; dashed line: v. The widths 
of the wave front and wave back (regions of fast change in u 
where (u, v) is not on the / nullcline) are of order y/e. The 
full return to equilibrium is not shown. 



for intracellular calcium dynamics; since the excitability 
here arises through the opening and closing of a small 
number of ion channels (allowing/preventing calcium ef- 
flux from calcium stores 0), fluctuations are inherently 
non-negligible. Indeed, experiments show direct evidence 
of noise effects in the form of abortive waves and spon- 
taneous wave nucleation[^|^]. 

In this paper, we study the process of spontaneous 
wave nucleation for a Id generic excitable system mod- 
eled by the Fitzhugh-Nagumo equations 
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Here r] v and r\ u are small independent white noise terms 
modeling fluctuation effects with covariance equal to: 

< r} i {x,t)r] j {x' ,t') >= f35 i:j 5{x - x')5(t - t') (3) 

At j3 = and for positive values of 7, this system is 
excitable with a single stable equilibrium point uq — 
— (l + y), vq = 3uq — Uq. As already mentioned, a wave of 
excitation can propagate through the system (see figQ); a 
counter-propagating pair of such waves will be generated 
if a local perturbation above a threshold value is applied. 
For negative values of 7, the system becomes oscillatory. 
This model is not meant to be a realistic approximation 
for any specific physical or biological process; instead we 
use the model to understand the generic features of wave 
nucleation due to noise. 

Specifically, at finite j3 the noise will allow the birth 
of pairs of counter-propagating wave at a rate that will 
depend on the amplitude of the noise. That rate can 
be determined in the case of relatively high noise using 
direct numerical simulations. However, in the case of 
low noise, such a method becomes computationally pro- 
hibitive. Here we present a computation of the transi- 
tion rate using a most probable escape path (MPEP) ap- 
proach that is based on solving the Fokker- Planck equa- 
tion 0, 0| . This method has been successfully applied 
to dynamical systems ^2 EH an d to a few cases of spa- 
tially extended systems, namely the transition from creep 
to fracture an d magnetic domain reversal To de- 
rive the equation for the MPEP, we use the fact that the 
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FIG. 2: Space time contour plot of the u (right) and «(left) 
field during nucleation with T=20; parameter values are e — 
0.1 and 7 = .2. Time increases upward and the white area 
corresponds to the region outside of [x m in(t) = max(xf(t) — 
x o ff,0),x max (t) = max(L + Xf(t) - x o ff,0)], [0,T]. The 
contours are at -2, -1, 0, and 1 for the u field and at vo — 0.2, 
Vo — 0.05, vo + 0.05, vo + 1, v a + 2 and v a + 3 for the v field. 
The lighter the surface, the higher the values of the fields are; 
thus, the light region of the u contour plot corresponds to 
the excited region. After the nucleation event, points near 
the center reach the maximum of the excited branch of the / 
nullcline (see fig. , return to the quiescent state and thereby 
give birth to the two wave-backs. 



solution of the Fokker-Planck equation with given initial 
and final field configurations for time interval [0, T] can 
be written in terms of the path integral 



P(T) = Jv(u,v)cxp (-^ J J dtdx 



S(x,t) (4) 



with the "action density" S given by the sum of the 
squared deviations of the time derivatives of the fields 
from their deterministic values: S = 5% + 6% 



S u = d t u — ('iu — u 3 — v)/e- 
5 V = d t v - u - (1 +7) 
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The functional integral is taken over all paths that be- 
gin at t — in equilibrium and end with a given final 
counter-propagating wave state (v,f, Vf) at t = T. Since 
we take j3 close to 0, the r.h.s. of eq. (@J is dominated by 
the path (called most probable escape path or MPEP) 
that maximizes the integrand over all paths; thus, the 
transition rate is found to be proportional to: 



exp(-£//3) 



(7) 



where E is the minimum over all paths of J dx J dt S 

Therefore, in order to compute the transition rate be- 
tween the rest state and a pair of counter propagat- 
ing waves, one has only to compute the minimum of 
the quantity in eq. Q). This minimum can be ex- 
pressed using a variational principle as the solution of 
a PDE 16]; however, using that approach to finding the 
actual MPEP between two different states involves the 
use of a shooting method with numerous parameters, 
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FIG. 3: (a) and (b) v profiles during a typical nucleation 
event, (a) v profiles at times 25, 25.5, 26, 26.5, 27, 27.5 
and 28. During that stage, the minimal value of the v field 
decreases, (b) v profiles at times 28(dashed), 28.25, 28.5, 
28.75 and 29 (solid lines). During this stage, the minimal 
value of v increases, (c): u profiles at the same times. The 
profiles corresponding to the (a),(b) curves are represented 
with a solid, dashed line, (d): The dashed line represents the 
/ nullcline. The thick solid line is the trajectory in the phase 
plane (u(x — 0,t),v(x — 0,t)) of the center of nucleation. 
The point where maximum value of J dxS^ + SI is reached 
is indicated by a thick arrow, whereas the point where this 
deviation is below 1% of this maximal value is indicated by the 
thin arrow. Beyond this point the dynamics are mainly driven 
by the deterministic equations. The thin solid lines represent 
spatial snapshots of the MPEP (u(x,t = Ti), v(x,t = T t )) 
profile at regularly spaced points in time (Ti — Tj_i = 0.5). 



which turns out to be numerically quite difficult. Instead, 
we used an alternative method based on discretizing the 
above path integral on a space-time grid and directly 
using a auasi-Newton[T^| method to find the minimum. 
One difficulty with this approach is that the Fitzhugh 
Nagumo model has a slow recovery time compared to the 
timescale associated with a pulse (width of pulse/speed). 
This then necessitates having a very large spatial do- 
main, if one attempts to completely encompass the re- 
gion over which the nucleated wave configuration differs 
from the quiescent fixed point. To get around this dif- 
ficulty, we use a grid moving with the pulse. That is, 
the grid moves in the same direction as the pulse in 
order to keep the wave front at a fixed distance from 
the boundary. Thus, the boundaries of the domain used 
to compute the MPEP path are no longer [0, L], [0, T], 
but [x mm (t) = + max(a;/(t) - af //,Q), x max (t) = 
L + max(ir f (t) — x Q f / , 0)] , [0, T] , where x / (t) is the posi- 
tion of the wave front defined as the point where u goes 
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FIG. 4: Parameters values are e = 0.1 and 7 = .2. 
Space time contour plot of the action density + 5% with 
its maximal value normalized to one. The contours are at 
1/96, 1/48, 1/24, 1/12, 1/6, 1/3 and 1/1.5. (b) Left: spa- 
tial distribution of the contribution to log 10 (_E) during the 
MPEP; Right: temporal distribution of the contribution to 
log 10 (_B) during the MPEP. One can see that the contribu- 
tion to £ as a function of time increases exponentially and 
then decreases rapidly to a very low level that is close to nu- 
merical noise 



above and x f / is an arbitrary value lower than L and 
significantly bigger than the wave front width. The fol- 
lowing boundary conditions are applied: 



at X — X ma x 
at X — Xjxiin 



(t): u = u a , v = v 

d x u(0) = if x min {t) = 
S u (0) = if (t)^o 
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The condition S u = implies that the recovery past x m i„ 
is purely deterministic and hence does not contribute to 
E. A check on our procedure is afforded by the fact that 
as long as the distance the wave had traveled (in the final 
state) is significantly bigger than the region over which 
the wave initiates, the results obtained are independent 
both of the distance the wave had traveled and of the 
width of the space window used (L). We fix T to be 
large enough that the deviation from deterministic dy- 
namics for very small time is negligible. Once T is fixed 
and the specific final wave state chosen, there is no time- 
translation invariance in the MPEP. 

We now describe the results obtained using this 
methodology. In the small e limit (e < 0.1), the shape 
of the MPEP is quantitatively independent of e and 7 
and even for high values of e, there is no significant qual- 
itative difference. In fig. [21 we present such a typical 
escape path. As shown on the contour plots, wave nu- 
cleation is found to be a very localized event. Essen- 
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FIG. 5: Left: circles log 10 e dependence of log 10 E using a log 
scale, solid line: y = | — 0.8. For small e, E behaves like 
yft. Right: circles 7 dependence of yfE for e = 0.1. Solid, 
y = 1.2 2: E behaves like 7 2 for a wide range of values of 7. 
The scaling breaks down for higher values of 7, close to the 
propagation threshold (for this value of e, the threshold is at 
7 = 0.5525 ±0.0005). 



tially, the noise acts to create a local dip in value of the 
v field which is then followed by a large positive excur- 
sion for the u field as it goes into the excited phase. To 
describe this mechanism more fully, we present in fig. [3] 
the phase-plane trajectory at the center of nucleation as 
well as several snapshots of the spatial form of the fields 
during the nucleation process. One can see then that 
the escape path consists of the center of nucleation being 
driven by noise below the minimum of the / nullcline; this 
then quickly drives the u field positive and leads after, v- 
field driven relaxation, to the pair of counter-propagating 
pulses. There is a significant fluctuation contribution to 
the nucleation event in a small region around the nucle- 
ation point(see fig. 0}. Results using different values of 
7 and e show that the width of that small region is pro- 
portional to the width of a front, that is yft (see fig. EJ. 
Note that the other simple possibility, that of nucle ating 
an excited region for the u field at a fixed value of t>|18|. 
is not observed. 

In accord with the shape of the trajectory, our calcu- 
lations show that the main contribution to E during the 
MPEP come from the J J d%, at least for small values of 
e. Thus for 7 = 0.2 the value of the ratio / / 5%/ f f 8% 
never went below 30 for e = 0.1, and went up to 1000 for 
the minimal value of e used (e = 0.001). For higher val- 
ues of e, the ratio was significantly lower (3 for e = 0.8, 
close to the limit of propagation for this value of 7 ) . Fur- 
thermore for small e, E scales like yfi (see fig. 01 • For 
higher values of e, this simple scaling is no longer valid. 
The fact that for small values of e, E scales like y/e and 
that a wave is therefore much easier to nucleate can be 
explained by a simple argument. As already mentioned, 
the phase-plane trajectory of the MPEP is mainly inde- 
pendent of the value of e. The only dependence is then 
due to the spatial scale appearing in the integral which 
is yfe, giving rise to the observed result 

We now describe results obtained when varying 7 with 
e held constant. This means that we consider the nucle- 
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ation rates for differing excitabilities but with the same 
time-scale ratio between the u and v field dynamics. For 
a wide range of values of 7, we find that E scaled like 
7 2 (see fig. EJ). The more excitable the system, the more 
likely it is to nucleate a wave. This scaling is that same as 
one would obtain analytically in the much simpler zero- 
dimensional version of this problem. Here the analog of 
wave nucleation is the noise-induced creation of a large 
transient excursion away from the fixed point. If e is 
small, the escape path will follow the / nullcline down to 
its minimum and then follow the noiseless dynamics to 
reach the other stable branch. In such a situation, one 
can compute the transition rate analytically in a piece- 
wise linear version of the FH model, 



-\{v + u) ifw<0, 
\(l-v-u) if m>0, 



•7. 
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After elimination of u, it can easily be shown that the 
variational equation for the MPEP has the form 



a v 

IF 



v + 7 



(11) 



with v(0) = vq = —7 and v(T) = 0. A straightforward 
calculation shows that the action is equal to 



7 



2(1 - exp(-2T)) 2 



(12) 



and that its minimum, reached for T — 00, is equal to 
7 2 /2. This result shows that the 7 2 scaling of E found 
in our computations can be interpreted as being due to 
7 equaling the "distance" of the stable equilibrium point 
from the border of its basin of attraction. Putting it 
all together, our data yield a log transition rate propor- 
tional to — 7 Y e . This result could be tested experimen- 
tally, perhaps by adding illumination noise to the light- 
sensitive BZ reaction. 

It is worth mentioning that there are other poten- 
tial applications of the MPEP approach to nucleation 
in spatially-extended non-equilibrium systems. One ex- 
ample concerns the thermal generations of localized 
patches of traveling rolls in electro-convection, as studied 
recently [Tij. Also, the method used here is not limited 
to white noise. A simple generalization of the derivation 
allows for the incorporation of multiplicative noise via 



dividing the (5 2 and 51 terms in the action density by the 
corresponding (possibly field dependent) variances of the 
noises added to the u and v equations respectively. This 
will clearly be necessary for the study of realistic calcium 
models. Finally, there is a similarity between the MPEP 
method and what must be done to consider quantum tun- 
neling in spatially extended systems [2^, where one also 
must find the entire space-time path tunneling trajectory 
in order to find a leading estimate of the rate. 
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